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A computer aided high temperature expansion of the mag- 
netic susceptibility and the magnetic specific heat is presented 
and demonstrated for frustrated and unfrustrated spin chains. 
The results are analytic in nature since the calculations are 
performed in the integer domain. They are provided in the 
form of polynomials allowing quick and easy fits. Various 
representations of the results are discussed. Combining high 
temperature expansion coefficients and dispersion data yields 
very good agreement already in low order of the expansion 
which makes this approach very promising for the application 
to other problems, for instance in higher dimensions. 
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I. INTRODUCTION 

Spin systems are among the most investigated systems 
in solid state physics. They represent problems with high 
correlation since the spin algebra does not have the sim- 
plicity of the fermionic or the bosonic algebra. This can 
also be seen from the generic derivation of antiferromag- 
netic spin models from a half-filled Hubbard model in the 
limit of large interaction U oo. Hence even the calcula- 
tion of simple properties like the magnetic susceptibility 
X or the magnetic specific heat C is not straightforward. 
Experimentally, however, susceptibility and specific heat 
are the first quantities used to characterise a compound. 
So quantitative theoretical predictions are very impor- 
tant to pinpoint the appropriate model. 

Quantum Monte Carlo methods underwent consider- 
able progress in the last years so that the calculation 
of X Slid C for unfrustrated systems has become feasi- 
ble to very high accuracy. The treatment of frustrated 
spin systems, however, is still a difficult task since in 
the standard Ising basis the sign problem occurs. This 
leads to the phenomenon that considerable cancellations 
and the concomitant loss of statistical accuracy occur in 
the quantum Monte Carlo computation, see e.g. Ref. [1]. 
In the case of strong frustration one still has to resort 
to complete diagonalisation which restrict the accessible 
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system sizes very much, see e.g. Ref. [2]. For (quasi) one 
dimensional systems like chains or ladders finite temper- 
ature density-matrix renormalisation provides a reliable 
means to calculate susceptibilities and specific heats, see 
e.g. Ref. [3-6]. The caveat, however, remains that the nu- 
merical methods require a new programme run for each 
set of parameters. So fitting becomes a tedious task as 
soon as more than one parameter is involved. 

The objective of the present article is to introduce a 
computer aided expansion in the inverse temperature for 
the quantities x ^-^d C. This method is a variant of the 
"linked cluster" approach [7, 8] . In the form implemented 
here it provides the results as polynomials in the rele- 
vant energy ratios. Thereby extremely fast and conve- 
nient fit procedures become possible. Frustrating terms 
do not pose more problems than any other additional 
couplings. For the sake of simplicity we will demonstrate 
our approach for one dimensional systems, i.e. chains. 
The subsequent choice of an appropriate representation 
of the results is also very important. The inclusion of 
low temperature information enlarges the range of valid- 
ity considerably. 

In the next two sections the method is explained in de- 
tail and contrasted to the conventional linked cluster ap- 
proach. In Sect. IV the results are given and represented 
in various ways in order to obtain the best description. A 
particular representation based on additional dispersion 
data is explained in Sect. V. The findings are summarised 
in the concluding Sect. VI. 



II. METHODS 

In the present work two methods are used to expand 
the physical quantities. The first one is the linked cluster 
method [7, 8] and the second one a method which will be 
called moment algorithm henceforth. To be explicit, the 
spin- 



■i Heisenberg chain 
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H = Y,{JSA+i + aJSA+2 - hSf) 



(1) 
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with nearest and next nearest neighbour interaction is in- 
vestigated with h representing the magnetic field in units 
of The ratio of nearest and next nearest neighbour 
exchange coupling is given by a. 

The main idea of the linked cluster method is to re- 
strict calculations to finite systems to obtain results in 
the thermodynamic limit. The method can be applied to 
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systems (clusters) which arc described by a sum of local 
hamiltonians like the Heiscnberg hamiltonian in (1). An 
expansion of a quantity in powers of such a hamiltonian 
results in the computation of the contributions of clus- 
ters of various sizes. To obtain the contribution of a finite 
cluster to the quantity of the infinite system all contribu- 
tions of subclustcrs have to be subtracted with suitable 
multiplicity. Only connected ("linked") clusters provide 
nonvanishing contributions. The non-conncctcd clusters 
cancel out due to the normalisation of the expectation 
value (see for instance Eq. (2)). This is the main result 
of the linked cluster theorem. 

The numerical approach of the moment algorithm 
makes use of the result of the linked cluster theorem. 
Here as well the physical quantities are evaluated in the 
thermodynamic limit by means of finite systems. Let us 
consider the Heiscnberg hamiltonian with nearest neigh- 
bour interaction acting on a chain and let us expand a 
physical quantity in powers of f3J = J /T, i.e. in pow- 
ers of this hamiltonian. The largest connected cluster at 
a given order n of the expansion contains {n + 1) sites. 
Hence all clusters of this size must be embedded com- 
pletely in the finite system to obtain valid results in the 
thermodynamic limit. 

For concreteness, we compute the magnetic suscepti- 
bility per site at vanishing magnetic field 



X{T) 



(3 tT{M^e-P^) 



(2) 



Denominator and numerator are computed separately by 
expanding the corresponding exponential functions. The 
resulting rational function is again expanded around P = 
to obtain a polynomial in the inverse temperature f3. 

On this stage the comparison of the moment algorithm 
to the linked cluster approach shows one advantage and 
one disadvantage. The advantage is that it is not neces- 
sary to determine and to classify all contributing clusters 
explicitly. This task may not be underestimated in view 
of the lack of efficient algorithms comparing graphs. This 
point matters in particular for complicated lattices with 
different types of bonds. In this respect, the moment al- 
gorithm is simpler than the linked cluster approach. The 
disadvantage is that the finite systems which have to be 
dealt with are fairly large, in particular for elevated or- 
ders in (3 and higher dimensions. In this respect, the 
linked cluster approach is better suited for higher dimen- 
sional problems. 

The disadvantage mentioned is less troublesome if an 
efficient way to compute traces of powers of the hamil- 
tonian is available. To this end we present an algorithm 
which computes such traces in a very fast way. "Fast" 
means that the necessary effort increases not exponen- 
tially with system size N but only in powers of N . The 
algorithm reduces the trace to an ordinary expectation 
value in a higher- dimensional Hilbert space. To this pur- 
pose the Hilbert space of the real system is doubled by 
introducing to each real site \ir) a "doubled" site 



Any operator A defined on the real Hilbert space acts on 
the tensor product of real and doubled Hilbert space in 
the canonical way. That is A becomes j4 1 acting as 
the identity on the doubled Hilbert space. 

Furthermore, we consider a state which is the (unnor- 
malised) product of singlets (or triplets with = 0) 
between the real and the doubled sites 



N 



\S) = lli\1rU)-\lrU))\ 



(3) 



i=l 



The key observation is that the trace in the original, real 
Hilbert space is identical to the expectation value of \S) 
in the extended Hilbert space 



tr(^)|,cal = I 

double 



(4) 



To see the identity (4) one considers first a single site. 
Explicit calculation shows 

tr(A)|,eai = (T 1^1 T> + a l^li) (5a) 

= (T4d \A\ Trid)U + {irU \A\ irU)\d (5b) 

= {S\A\S)U (5c) 

where we used the subscripts r and d for 'real' and 'dou- 
bled' (extended) Hilbert space, respectively. The validity 
of (4) follows from (5) directly for all operators A which 
are products of local spin operators since the Hilbert 
space of many spins is just the tensor product of the 
Hilbert spaces of the individual spins 

trJjAil^ = JJtr^il^ (6a) 

i i 

= li^iilrid I - {irU \)Aii\ Trid) - | irU))\i (6b) 

i 

= {S\]lA,\S)U. (6c) 

i 

/,From the linearity of the expectation value and the trace 
follows then that (4) holds also for all operators A since 
they can be decomposed into sums of products of local 
spin operators. 

On the left hand side of (4) for A = for instance one 
has to compute in the Ising basis (spins up or down) 2^ 
contributions for terms {N number of sites, L number 
of bonds), i.e. one has to sum about L^2^ terms. On the 
right hand side of (4), however, one starts with a single 
state which will be excited in L ways by the application 
of H (g) 1 and requires to be de-excited in the same way 
so that only L terms contribute in the end. Thus one 
saves an exponential factor. An obvious by-product are 
the relations 



(7a) 
(7b) 
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which imply that for a given order to in /3 one needs to 
calculate only about to/2 applications of H to the singlet 
product state \S). This statement remains true for the 
numerator of susceptibilities as in (2) if the observable 
(here: M = ^t) commutes with H . This is the 

case for the uniform magnetisation M. Replacing \S) by 
MIS') thus makes the relations (7) also applicable to the 
numerator of (2). 

We implemented the actual calculations on computer. 
For not too high orders this can still be done with com- 
puter algebra programmes. But to obtain the highest 
orders it is necessary to write task-specific programmes. 
This was done by using the language C"'""'". Yet the de- 
pendencies on all coupling constants are included on the 
symbolic level, i.e. in polynomials of the coupling con- 
stants. So, once obtained, the results are available to 
everybody and they can be fitted to any experimental 
curve very quickly and easily. The results for the mo- 
ment algorithm (frustrated (order 10) and unfrustratcd 
spin chain (order 16)) and for the linked cluster algorithm 
(unfrustrated spin chain only, but order 24) are presented 
in the appendix. 

Before concluding this section we like to note that the 
trick to pass from a trace to an expectation value in a 
higher dimensional Hilbert space is not restricted to spin- 
i systems. By introducing for instance the generalised 
product state of singlets \S)g 

2S 

\S),=\[Y,^^^\{S-^)rA^-S)a) . (8) 

the method can be applied to arbitrary spin. In (8) ^d) 
stands for the state where the real spin has = j and 
the doubled spin = I. For a derivation one simply has 
to redo the calculation (5). 

Henceforth the angular brackets will denote the ex- 
pectation value with regard to \S) or the original trace, 
respectively, since their identity is established and so no 
further distinction necessary. 

The calculations in the present work were done to the 
highest order possible on the available work stations. For 
the unfrustrated chain the physical quantities were ex- 
panded to order iV on a finite system of size N sites. 
Therefore clusters in the A''th order are overcounted or 
missed. But it is possible to correct these effects in high- 
est order by an analytical argument which is presented 
in the subsequent section. 

III. FINITE SIZE CORRECTIONS 

The wrap-around effects of the numerator and denom- 
inator of Eq. (2) are investigated separately for a chain 
of length N with periodic boundary conditions. 

The denominator has the following kind of contribu- 



tions in the A'^th order 

N 

{H'') = {l[{SA+^)) + ... (9) 

i=l 

which are not realised in the thermodynamic system. 

Fixing the component of one of the spin vectors involved 
to for instance one sees that a non-vanishing contri- 
bution occurs only if all spin components are S*^. The 
overall value of the right hand side in (9) is 4~^. Since 
all permutations of the scalar products will occur as well 
if the left hand side of (9) is expanded and since all these 
permutations yield the same contribution the factor A^! 
has to be added. Since the choice of the spin compo- 
nent was arbitrary an additional factor 3 concludes the 
argument. Thus one has to subtract 3iV!4~^ to yield the 
thermodynamic result in the denominator. 

The corrections of the numerator of Eq. (2) are more 
complicated. They consist of three contributions. The 
first is similar to the one in the denominator 

N 

{{StfH'') = {{S^r UiSA^,) + ... (10) 

i=l 

and overcounts the numerator by 3iV!4~^~^. On the 
other hand, the thermodynamic contribution is neglected 
in this order of expansion 

{S!Sf,+,{Sj2){S2S3) . . . {SnSn+i)) (11) 

which represents the second correction. It takes the value 
27V!4~^~-'^. The factor A'^!4~^~-'^ arises for the same rea- 
sons as before. There is no factor 3 since the spin compo- 
nent is already fixed by the choice of the magnetisation 
direction. The geometric factor 2 arises instead because 
the sites 2 to A'' + 1 can be found to the right or to the 
left of the starting site 1. 

The third and last correction consists of clusters with 
triple occurrence of two sites. It is based on the identity 
for 5 = 1/2 

S^SyS' = i/8 (12) 

which holds also for all cyclic permutations of the spin 
components. For anticyclic permutations the left hand 
side of (12) acquires a minus sign. 

A wrap-around with sites occurring three times is pos- 
sible as soon as the magnetisation operators are taken to 
be at different sites: SfSj with j ^ 1. Then each of the 
sites 1 and j appear three times 

N 

{S!S]l[{SA+^) . (13) 

i=l 

Note that all permutations of the sequence of the scalar 
products appear. Hence in most cases the contributions 
cancel each other since cyclic and anticyclic permutations 
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of the spin components at site 1 or site j occur indepen- 
dently and with equal amplitude. Only if the site j is 
adjacent to 1, i.e. j = 2 or j = N (periodic boundary 
conditions), the cyclic and anticyclic permutations at site 
1 and j are correlated and a finite total effect remains. 
The relevant factors are (for j = 2) 



{StSliSNSi){SiS2){S2S3)...) 



(14) 



Among the 3! = 6 ways to arrange the three scalar 
products in (14) the two where (5152) is in the mid- 
dle yield 1/8^ (5^ 5f + 5^5|) while the other four yield 
-l/82(5^5f + 5^5|) so that -2/8'^{S%S^ + 5^5|) re- 
mains. Accounting for the multiplicity due to the ar- 
rangement of the other scalar products yields the combi- 
natorial factor N\/3\. A factor 2 comes from the possibil- 
ity to choose j = 2 or j = N. The overall third correction 
finally reads -8 • 4-^+^N\/3\. 

In summary, the total corrections to the results com- 
puted in the iVth order for a finite system of A*" sites with 
periodic boundary conditions for the numerator Nu and 
for the denominator De are 



■^^^corrected ^computed A^! * ^ I ^ 



AT -1-1 



AT 



(15a) 



Decorrected = Decomputed - A'! • 3 ^- J . (15b) 

After these considerations it is also straightforward to 
correct the wrap-around effects for the Heisenberg chain 
with next-ncarost neighbour interaction. In the A'th or- 
der on a finite system of 2A^ sites these effects occur only 
in the Nth order in a and in f3. In this order the sys- 
tem corresponds to a system of two independent chains 
with nearest neighbour interactions only, so that the cor- 
rections (15) apply to terms of A'th order in a. This 
concludes the discussion of the finite size corrections. 



IV. RESULTS AND REPRESENTATIONS 

In the appendix the series coefficients are presented 
for the magnetic susceptibility x ^-iid for the magnetic 
specific heat C per site. The specific heat is derived from 
the denominator of Eq. (2) which is the partition function 
of the system by 



C{T) = 




(16a) 
(16b) 



It is worth mentioning that due to the derivation in (16b) 
one order in /3 is lost. It is re-gained, however, by the 
subsequent derivation with respect to T. 

In particular, results for the unfrustrated chain are 
listed in Appendix 3. These are obtained by the linked 



cluster method and comprise orders as high as 24. Such 
high orders are obtained by an exact extrapolation of 
smaller clusters as explained below. They are in agree- 
ment with expansion results previously obtained [9, 10]. 
In Appendix 3 a the magnetic susceptibility coefficients 
are given, in Appendix 3 b the specific heat coefficients 
are given. 

For the loop-free chain the calculation of high tempera- 
ture series for spin models simplifies decisively compared 
to higher dimensions due to the simple structure of the 
contributing cluster. A finite open chain with n bonds 
will contribute to the specific heat only in order 
This is so because otherwise there will always be one site 
occuring an odd number of times in the spin products 
leading eventually to a vanishing trace. 

The susceptibility series shows a systematic pattern 
which can be exploited to obtain longer series. For the 
S = 1/2 Heisenberg chain, the nonvanishing contribution 
of order /3" of a finite open cluster with n bonds and n+1 
sites will be due to one spin product of the hamiltonian 
acting on each bond and two magnetisation operators 
at the ends of the chain. Its contribution is of the form 
agn! 4~"~-'^/3", where gq is independent of the cluster size 
(as before we consider the terms without the 1/n! factors 
from the exponential series). 

Similarly, the contribution in order will be of the 
form {bo -I- n&i)n! 4~"~^/3"+^. Again bo and bi do not 
depend on the length of the finite cluster. The term 
proportional to n is due to the n possibilities to attach 
one more term of the hamiltonian to any of the n bonds. 

In the same way, the general form of the higher order 
contributions can be determined with more and more co- 
efficients. With sufficiently large clusters the coefficients 
ao, (6o,6i), (co,ci,C2), ... can be obtained by solving a 
system of linear equations. Using clusters with up to 18 
bonds allowed to extend the susceptibility series to order 
24. 

The moment algorithm allowed us to compute results 
for the unfrustrated chain up to order 16 for the sus- 
ceptibility (Ala) and for the specific heat (Alb). Note 
that this is only two orders less than the maximum clus- 
ter which is actually computed in the linked cluster ap- 
proach. The results are in complete agreement with the 
linked cluster results. The only difference is in the inter- 
nal representation where doubles are used in the linked 
cluster programme and true fractions in the moment al- 
gorithm. The susceptibility as well as the specific heat 
of the frustrated chain is expanded up to order 10 (A 2 a 
and A 2 b). 

Having obtained the series for the various expansions 
we pass now to the discussion of suitable representations. 
The choice of an appropriate representation allows to 
gain the maximum of information from the bare series 
coefficients. We follow two main routes. One is the use 
of Fade approximants and continued fractions, respec- 
tively; the other is to incorporate additional information 
at low temperatures to improve the representations in 
the low temperature regime. 
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Exact results from Bethe ansatz calculations for the 
unfrustrated chain [11, 12] and numerical results from 
density-matrix renormalisation group (DMRG) calcula- 
tions for the frustrated chain [13] are used as benchmarks. 



A. Unfrustrated Chain 

Since the convergence of the plain series in /3 can be 
hindered by any pole the use of a Pade approximant de- 
scribing the quantity under study by a rational function 
is more stable than the plain series. 
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Figure 1. Comparison of various depths of continued frac- 
tion representations of the susceptibility x for the unfrus- 
trated chain. 



The polynomial in (3 of the physical quantity under 
consideration is represented as a continued fraction, say 
with depth 2N , 



Cl + 



(17) 



C2 



C2N 



which is equivalent to the [A'^, A'^] Pade approximant. An 
odd depth of 2Ar-|-l is equivalent to the [iV -|- 1, A''] Pade 
approximant. Increasing the degree of cither the numer- 
ator or the denominator polynomial at the expense of the 
other does not improve the results. The advantage of the 
continued fraction representation is that the coefficients 
Ci remain constant on increasing depth. For instance ci 
is equal to 4 and C2 is equal to 1/2 due to the Curie and 
due to the Curie- Weiss law, respectively. 

In Fig. 1 various depths of continued fraction represen- 
tations of the unfrustrated susceptibility are shown. The 



comparison with the exact result shows that the agree- 
ment improves on increasing depth as expected. But 
the improvement of the representations is relatively small 
for higher depths. Excellent agreement can be achieved 
down to T « J/4 if coefficients up to order 24 are used. 

In order to extend the region of satisfying agreement 
to lower temperatures information about the low tem- 
perature regime can be incorporated. To this end, the 
continued fraction depth is incremented by adding a new 
constant. This constant is not determined from the high 
temperature expansion. But it is determined such that 
the desired additional property is fulfilled. 

To be precise, we include the value of the susceptibility 
at zero temperature. For the unfrustrated chain it can 
be expressed as [11, 14] 



X(0) = 



1 1 

27r vs 



(18) 



with the spin wave velocity vs = ^ [15]. Eq. (18) implies 
that the central charge c is one. Assuming that the cen- 
tral charge does not vary on switching on the frustration 
we will use (18) there, too. For the change of the spin 
wave velocity is accounted by [16] by 



vs 



-(1 - 1.12a) for < a < Qfc 



(19) 



In the gapped regime for a > ac « 0.241167 [17] the 
susceptibility vanishes exponentially at T = 0. 

The relevant gap, however, is not the spectroscopic gap 
Aoi between the S = ground state and S = I excited 
states but half of this value Aoi/2. This is so since the 
elementary excitations of strongly frustrated spin chains 
are asymptotically free massive S = 1/2 spinous, see for 
instance Refs. [18-20]. 

The low temperature behaviour of the specific heat [11] 
is given by 



C(T « 0) 



TT 1 

3 vs 



■T 



(20) 



with the same spin wave velocities vs as in the previous 

equations 

for a < ac- In the gapped regime for supercritical frus- 
tration, an exponential vanishing for low temperatures is 
to be expected. Prom (20) follows directly 



;^^(^ = o) = ?^ 

al 3 Vs 



(21) 



which can also be incorporated in the representations. A 
third piece information is obtained by 



s(oo) - s(0) = 



cm 



dr = ln2 



(22) 



This piece of information, however, is more difficult to 
build-in since it involves an integration over the contin- 
ued fraction. Moreover, it turns out that its effect is not 
sizable. Thus it is not considered any further. 
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Figure 2. Various representations of the susceptibility x 
for the unfrustrated chain. The plain series is expanded up 
to order 24. The inset shows a zoom of the [12,ll]-Paxle ap- 
proxirnant and of the [12,12]-Pade approximant with T = 
information, see Eq. (18). 



Fig. 2 shows the various representations for the sus- 
ceptibility X of the unfrustrated chain. The approxi- 
mate results agree very well with the exact ones down 
to T/ J pa 0.2 for the Pade approximant with T = in- 
formation. Without the aid of the exact result one is 
also able to determine the quality of the representation 
by comparison of results in highest order with those in 
lower orders. 



estimate for the radius of convergence of the plain se- 
ries. From the values given in Fig. 3 one can deduce that 
this value is fairly constant at about /3max ~ 1.9. This 
implies that the plain series will always diverge below 
about T K. 0.53 irrespectively of the order of the series, 
cf. Figs. 1, 2. This is a good illustration of the utility 
of Fade representations. They are not blocked by the 
occurrence of poles. So they are able to represent more 
complicated functional dependencies. 




Figure 4. Various representations of the specific heat C for 
the unfrustrated chain. The plain series is expanded up to 
order 24. The inset shows a zoom of the [ll,ll]-Pade ap- 
proximant and of the [12,12]-Pade approximant with T = 
information, see Eq. (20) and (21). 
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O [8,8] Pade, smallest modulus: 1.94 
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Figure 3. Singularities of various Pade representations in fi 
for X and the smallest moduli of their /3 values. Singularities 
in the left half-plane are not shown. 



It is instructive to look at the poles of the Pade rep- 
resentations closest to the origin. Their modulus is an 



The low temperature behaviour of the specific heat is 
less complex than the one of the susceptibility. Fig. 4 
shows a better agreement with the exact result, especially 
in the low temperature regime. Here the plain series is 
expanded up to order 24 in /3. With the two pieces of 
low temperature information (20,21) the representation 
agrees very well with the exact result down to T/ J « 0.1. 

At this stage a consideration of the range of valid- 
ity that one could expect is in order. In fact, wc argue 
that one should have expected an even better description 
based on 1/T results. 

Calculating up to order n means that the physics on 
a length scale n (lattice constant set to unity) is taken 
into account since this is the size of the maximum cluster 
treated properly. So one is led to the estimate 



vs 



(23) 



where the energy scale 27rrmin results from the discreti- 
sation of the Matsubara frequencies which serves here as 
infrared cutoff. From (23) follows for the unfrustrated 
chain Tmin 1/ (4n) for the temperature down to which 
the large T information should be capable to describe the 
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physics properly. It is obvious that the vahdity stops ac- 
tually at much higher temperatures. For this reason we 
presume that the representation by a Pade approximant 
is not yet the optimum. 

B. Frustrated Chain 

Motivated by the inorganic spin-Peierls system 
CuGeOs [21,22] the results for the frustrated chain are 
presented with a fixed a- value of 0.35 [23-25]. This value 
is chosen since it allows a good description of the suscep- 
tibility data. At low temperatures there is evidence that 
the frustration is lower [26]. 

Comparisons to benchmark calculations were also per- 
formed at the critical frustration ac- Compared to the 
higher orders reached in the unfrustrated case, the re- 
sults for the frustrated chain should agree well only for 
higher values of T/J. We will see, however, that this is 
not the case. 
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Figure 5. Various representations of tlie susceptibility x 
for the frustrated chain with a = 0.35. The plain series is 
expanded up to order 10. The inset shows a zoom of the 
[5,4]-Pade approximant and of the [5,5]-Pade approximant 
with T = information for the gapped regime. 



we had to include much higher orders to achieve similar 
agreement down to T/ J w 0.25. The Figs. 5 and 6 depict 
data for a specific value of frustration. But the raw data 
as given in Appendix A2 allow the calculation of suscep- 
tibilities and specific heats for any value of frustration. 

Three possible sources for the improvement of the 1/T 
description by frustration are conceivable. One is the ap- 
pearance of a gap due to frustration. But for a = ac we 
found qualitatively the same behaviour so that this ex- 
planation can be excluded. A second idea concerns the 
dominance of logarithmic corrections. Since our ansatze 
are not fit to represent these corrections the agreement 
must deteriorate once logarithmic corrections become im- 
portant on lowering the temperature. If this mechanism 
were the dominant one one should expect a significantly 
improved agreement at the critical frustration. The ac- 
tual comparison (not shown), however, does not display 
a significantly improved agreement. So the logarithmic 
corrections seem to be not the main problem of a correct 
representation [27]. 

The third possible explanation is a reduction of the 
spin wave velocity or, put differently, of the whole dis- 
persion. Analytically, it is known in leading order of an 
expansion around the dimer limit that frustration low- 
ers the mobility of the excitation [28] . Numerical results 
show the same, see Eq. (19). Indeed, the positions of 
the maxima and the lower bound of the range of valid- 
ity scale roughly like the spin wave velocity as given by 
Eq. (19). Hence, our results indicate that the estimate 
(23) is valid to the extent that it establishes a propor- 
tionality Tniin vs/n. 

Summarising this section we state that the Pade rep- 
resentations with low temperature information incorpo- 
rated show very good agreement down to rather low val- 
ues of T/J. In particular the maxima of the physical 
quantities susceptibility and specific heat are sufficiently 
well described. With the full dependence of the model 
parameter one has a powerful tool to fit the parameters 
to experimental data in a very fast and convenient way. 



Fig. 5 shows the susceptibility compared to DMRG 
calculations. The best representation with T = infor- 
mation is in very good agreement down to T/J 0.25. 

Since the frustration is supercritical 0.35 > ttc the 
T = information consists in fixing x(r = 0) = due 
to the exponential vanishing. The region of satisfying 
agreement coincides very well with the one of the repre- 
sentation of the specific heat C with T = information 
in Fig. 6. For the specific heat in the supercritical regime 
the derivative dC/dT is set to zero, too. 

Obviously, frustration is favourable for the range of 
applicability of the 1/T expansion. Without frustration 
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U 0.2 




Figure 6. Various representations of the specific fieat C for 
a frustrated chain with a = 0.35. The plain series is expanded 
up to order 10. The inset shows a zoom of the [4,4]-Pade 
approxirnant and of the [5,5]-Padc approximant with T = 
information, see Eqs. (20) and (21) 



V. REPRESENTATION WITH DISPERSION 

DATA 

In this section a different kind of representation is il- 
lustrated. It is also based on the idea to incorporate 
low temperature information in a high temperature ex- 
pansion. As in the previous section, the approach is 
motivated by approximations for the susceptibility of a 
dimerised and frustrated S = ^ chain in Refs. [29,30]. 
There the authors approximate the magnetic susceptibil- 
ity by using an exclusion statistics appropriate for excited 
triplets in the dimer model. To first approximation these 
excitations are often treated as free bosons. Yet it is ob- 
vious that they are hard-core bosons since there cannot 
be more than one on each dimer. This basic fact was 
included in Ref. [29]. In Ref. [30] the interaction beyond 
the hard-core exclusion was incorporated on a mean-field 
level. 

Defining the partition function 



271" J-^ 

for a single excitation one obtains the hard-core 



Xo=/3 



z{(3) 



1 + 3z(/3) 



(24) 



(25) 



for the hard-core exclusion statistics. On the mean-field 
level one obtains 



X 



Xo 



1 + -'^effXO 



(26) 



where Jcfr can be determined either in the limit of strong 
dimerisation or in such a way that the Curie- Weiss con- 
stant is correct. Both methods do not differ much [30]. 
The approach (26) is very successful in describing triplets 
with small dispersion [31] as in SrCu2 (603)2 [32]. 

Formulae (25,26) suggest to represent Tx essentially as 
function of z{(3). In the lowest order the representation 
should reproduce Eq. (25). Furthermore, it should allow 
to incorporate the information of the high temperature 
expansion in a natural way and the ansatz should be as 
simple as possible. Our choice is 



coz{f3) 



1 + 



civ{l3) 



(27) 



1 + 



1 + C3V{I3)' 



with 



v{(3) = 1 - z{(3) . 



(28) 



The variable u(/3) is chosen such that w(/3) cx /3 for 
/3 ^ so that the coefficients in (27) can be determined 
straightforwardly from the high temperature expansion. 
Of course, the choice (27) is just one of many possible 
choices so that a certain degree of arbitrariness remains. 
We tried also other choices, for instance an ansatz ex- 
tending (26) where xo is taken as variable instead of v. 
Our observation is that the particular choice does not 
matter much so that we present here the easiest ansatz 
we could think of. 

The low temperature information incorporated in (27) 
is in the dispersion relation uj{k). In order to demonstrate 
the approach we apply it to the unfrustrated chain where 
we can rely on exact results for the dispersion [15] 



uj{k) = — sin(fc) 



(29) 



We are aware that the unfrustrated case is not partic- 
ularly suited for the approach (27) as motivated above. 
The elementary excitations are 5=1/2 spinous [33], not 
magnons. In this respect, we are choosing a difficult test 
case for which we will show that the approach works very 
well. On the other hand, it is known that the main weight 
of the dynamic structure factor [14, 34] is located close to 
the lower boundary given by (29) so that the use of (29) 
as "magnon dispersion" is justifiable. 
Evaluation of the integral (24) yields 



z{(3) = lo ( i/?^] - Lo (^(37T 



(30) 



with the modified Bessel function of the first kind Ii, and 
the modified Struve function as defined in Ref. [35]. 

By construction, the ansatz (27) is able to fulfill the 
high temperature limit /? — > where Tx — > 1/4. It does 
so if Co is set to 1/4. It is a very favourable feature 
that the opposite limit of vanishing temperature /? ^ 00 
where Tx — > T/n^ [36, 37] can also be reproduced. Using 
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z{0) = 2/7r 



7r/2 



-/37r/2sin(/c)* 



(31a) 



poo 

2/77 / e-'^^/^fc^ ^ 4/(^7r2) (31b) 

JO 



one easily sees that the correct T — > limit is obtained 
if 



1 = 4co/(l + ci/(l + C2/(l + cs/ • • • ))) (32) 



holds. 



0.15 



?< 0.12 




Figure 7. Various orders of continued fraction representa- 
tions of Eq. (27) for the magnetic susceptibility of a Heisen- 

berg chain in comparisou to the exact result obtained by 
Bethe ansatz. The arrow indicates the exact result at T = 0. 



Fig. 7 shows various continued fraction representations 
of the magnetic susceptibility of a Heisenberg chain. Al- 
ready low orders of the representation show good agree- 
ment with the exact result down to low T/J. Even the 
sixth order representation describes position and height 
of the maximum fairly well. For order 8 and above, the 
maximum is perfectly described. Even the value of the 
zero temperature susceptibility is very close to its exact 
value which could not be expected since the validity of 
(32) is not built-in. 

The 12th order representation fits very well down to 
T/J « 0.3, which is almost the result of the [12, 12]-Pade 
approximant in Fig. 2 obtained from a much costlier se- 
ries up to order 24. Orders above 12 are difficult to im- 
plement with the dispersion information since the deter- 
mination of the constants Cj becomes very tedious. Par- 
allely, the improvement obtained becomes smaller and 
smaller. 

Looking at Fig. 7 closely it can be concluded that only 
the logarithmic terms in the susceptibility [12, 38] spoil 
the agreement at low temperatures T/J < 0.3. If these 
terms are explicitly incorporated the agreement in the 



whole temperature range can be conveniently described 
[6]. 

It is, however, not our aim to provide a fit to a result 
which is known analytical. For this reason we do not 
follow the route to incorporate the logarithmic terms into 
the ansatz (27). By the results depicted in Fig. 7 we have 
demonstrated that the inclusion of T = information in 
an ansatz of high temperature expansion improves the 
range of validity considerably. In particular, already a 
small number of high temperature coeSicients allows a 
satisfyingly accurate description of the overall form of 
the physical quantity under study. To corroborate this 
conclusion we present in Fig. 8 the analogous result for 
the specific heat. It is based on the ansatz 



do (z"-3(zOV(l + 3z)) 



(33) 



1 + 



d2v{P) 



where z' and z" stand for the first and the second deriva- 
tive of z with respect to /3, respectively. The ansatz (33) 
is motivated by the result 



{z" -3{z')y{l + 3z)) 
l + 3z 



(34) 



derived from the free energy including exclusion statis- 
tics [29] . As for the susceptibility the agreement between 
approximate ansatz and exact results is good even for 
low depths of the continued fraction (33). Note, how- 
ever, the spurious pole at about T = 0.3J occuring in 
the representation of depth 6. This phenomenon cannot 
be excluded in Fade representations so that one should 
always consider various depths in order to judge which 
features are meaningful. The 12th order result agrees ex- 
cellently with the exact result which shows the efiiciency 
of the ansatz (33). 



0.40 



0.30 



U 0.20 



0.10 



0.00 
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Figure 8. Various orders of continued fraction representa- 
tions of Eq. (33) for the specific heat of a Heisenberg chain in 
comparison to the exact result obtained by Bethe ansatz. 

Against the ansatze (27,33) one may object that the 
dispersion u!{k) will not be available in general. Of 
course, exact results for the dispersion are as rare as exact 
results for susceptibilities or specific heats. But there are 
a number of good approximate methods which provide 
Lu{k) at zero temperature, for instance perturbative ex- 
pansions [39] . Their results can be taken to refine and to 
supplement the high temperature expansions. Or it can 
be reasonable to use the experimentally determined re- 
sults for Lu{k) in order to understand the thermodynamic 
quantities. 

The fact that already low orders of a high temperature 
expansion can be sufficient to provide good estimates for 
x{T) and C(T) is especially important for higher dimen- 
sional models in d = 2 or d = 3. In these cases higher 
orders cannot be obtained due to the quickly rising num- 
ber of sites or number of clusters to be considered. 



VI. SUMMARY 

The main objective of the present paper is to illus- 
trate a general way to obtain very fast and convenient 
analytical formulae for standard thermodynamical quan- 
tities such as the magnetic susceptibility and the specific 
heat. So the results should be viewed mainly as effective 
tools for quick data analysis. The route to such formulae 
comprises two steps. The first one is a high temperature 
expansion in P performed symbolically on computers pro- 
viding the coefficients in analytic form. The second one 
is the use of an optimised representation of the results. 
In this second step the inclusion of additional informa- 
tion at zero temperature available from other sources is 
particularly useful. It is not our objective to compute by 
the methods presented unkown low temperature physics. 

For the sake of illustration we considered in the present 
article frustrated Heisenberg chains, i.e. chains oi S = 
1/2 spins with nearest and next-nearest neighbour cou- 
plings J and aJ, respectively. For these chains we pro- 
vided the high temperature coefficients for the magnetic 
susceptibility and the specific heat up to order 10 in the 
frustrated case and up to order 24 in the unfrustrated 
case. 

For topologically simple lattices the linked cluster ap- 
proach is the most efficient since it avoids to compute 
powers of the hamiltonian on unnecessarily large sys- 
tems. The actual expansion is done for subsystems, the 
so-called clusters. The price to pay is the necessary so- 
phisticated bookkeeping of the clusters. 

In order to be able to treat straightforwardly also more 
complicated lattices (or more complicated topologies of 
couplings such as frustrating couplings) we abandoned 
the calculation on subsystems in the moment algorithm. 



In order not to be overwhelmed by the quickly rising di- 
mension of the Hilbcrt space we identified the original 
trace with an expectation value on an extended Hilbert 
space. This trick reduces the number of terms in the 
application of the hamiltonian from L2^ to L where L 
is the number of bonds and the number of sites. In 
the moment algorithm the inclusion of frustration is not 
more complicated than the inclusion of any other addi- 
tional coupling. This is in contrast to quantum Monte 
Carlo approaches where frustration leads generically to 
the severe sign problem. 

For not too low temperatures a very good agreement 
could be achieved. Including further information on the 
low temperature the agreement can be improved further. 
In particular, the use of information on the zero temper- 
ature dispersion turned out to be very efficient. In this 
way, even relatively low orders allow a good description 
of the thermodynamic quantities under study. The nec- 
essary dispersion information can be taken from exact 
or approximate results. Taking experimental dispersion 
results allows to check the consistency of the model as- 
sumed. 

Due to the possibility to reach satisfactory results al- 
ready in low orders of the high temperature expansion 
the application to higher dimensional cases such as the 
strongly frustrated Shastry-Sutherland model [40,41] is 
possible and will be reported elsewhere. Work on gapped 
systems such as dimerised chains is in progress. 
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APPENDIX A: COEFFICIENTS 
1. Unfrustrated chain, moment algorithm 



Here the results for the unfrustrated chain are presented. The susceptibihty and the specific heat have been 
computed up to order 16. 



a. Susceptibility 



h. Specific heat 
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Scries coefficients a,i for the high temperature expan- 
sion of the magnetic susceptibility x = ^ Ylm ^n{PJ)"- 
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Scries coefficients a„ for the high temperature expan- 
sion of the magnetic specific heat C = J2n ^niPJ)"'- 



2. Frustrated chain, moment algorithm 

The coefficients for the results of frustrated chain are presented. The magnetic susceptibility and the magnetic 
specific heat are expanded up to order 10 in /3J. 



a. Susceptibility 
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Series coefficients a„ k for the high 



temperature expansion of the magnetic susceptibility of the frustrated chain 



12 



b. Specific heat 
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Series coefficients a„ k for the high temperature expansion of the magnetic specific heat of the frustrated chain 

3. Unfrustrated chain, linked cluster expansion 



a. Susceptibility 



n 


an 


n 


an 


n 


a„ 


n 


an 


n 


an 





1.0 


5 


-4032.0 


10 


-9565698560.0 


15 


-205019990184689664.0 


20 


-18366266410738921187573760.0 


1 


-4.0 


6 


-89376.0 


11 


-210597986304.0 


16 


-3169755454477500416.0 


21 


-40780317289246872850923520.0 


2 


0.0 


7 


163840.0 


12 


3486950684672.0 


17 


208763541109969256448.0 


22 


38668138493195891009425244160.0 


3 


64.0 


8 


26313984.0 


13 


203634731188224.0 


18 


8342101010835559022592.0 


23 


983734184997038611238624428032.0 


4 


400.0 


9 


191334400.0 


14 


-127324657152000.0 


19 


-175912858271144581529600.0 


24 


-75650797544886562610211717119286.7 



Scries coefficients for the hnked cluster expansion of the magnetic susceptibihty for the Heisenberg chain with % = 

^ an ( J \n 

4T Z^n (n+iy.yiTl " 



b. Specific heat 



n 


an 


n 


an 


n 


an 


n 


an 


n 


an 





0.0 


5 


-7200.0 


10 


-11558004480.0 


15 


-121944211136778240.0 


20 


96147483542540314214400.0 


1 


0.0 


6 


15120.0 


11 


199812856320.0 


16 


8791781390116945920.0 


21 


1279121513829538179364945920.0 


2 


6.0 


7 


1848672.0 


12 


10106191180800.0 


17 


310402124957945954304.0 


22 


27962069861743501862336200704.0 


3 


36.0 


8 


11426688.0 


13 


-19376365252608.0 


18 


-7225535925744106143744.0 


23 


-2398518627113966015427501883392.0 


4 


-360.0 


9 


-594846720.0 


14 


-9289795522775040.0 


19 


-643407197363813620776960.0 


24 


-129834725539335848980192847460554.1120 



Series coefficients for the hnked cluster expansion of the magnetic specific heat for the Heisenberg chain with C = 

Ear,. ( J \n 



13 



